/*************************************************************************************************************/
/* 06-05-2003: Ver 4; this is the first functional vertion of the fastICA program, related to ver 3,         */
/* the changes are:                                                                                          */
/*		    - variables that contain the calculations were renamed in order to agree with the        */
/*                    fastICA documentation available.                                                       */
/*                  - the variables mentioned above were changer from private to public in order to          */
/*                    provide data access to graphic modules.                                                */
/*                  - some limited access to graphic methods were tested successfuly with em_graph class.    */
/* 14-10-2004: Program changed name from InputIcaData to Ica preprocess. Began efforts to integrate ICA      */
/*                  to the hole portfolio program.                                                           */
/* 27-03-2005: Statistical and non escencial methods that can be reused where moved to the miscellanea class */
/* 28-03-2005: The name of the class waschanged again to Preprocess in order to show that the class has the  */
/*             capability to perform both: PCA & ICA preprocess.                                             */
/*************************************************************************************************************/
package preprocess;

import java.io.*;
import java.lang.Object.*;
import java.math.*;
import Jama.*;

import misellanea.Miscellanea;


public class Preprocess
{  
        int iColumn;
        int iRow;
        private double [][] dDataTrain = null;          // Input data train,
        
	private static int  iNComp;			// Number of components to be extracted.
        private static int  iMaxIter;                   // Max number of iterations
	private double [][] dCentMtx     = null;
	private double [][] dScaledMtx   = null;
        //private double [][] dDvals       = null;         // EigenValues (SVD)
        private double [][] dEigenValues = null;         // EigenValues 2
        private double [][] dKMtx        = null;
        private double [][] dWMtx        = null;
	private double [][] dX_Mtx       = null;        // Pre-processed MAtrix.
	private double [][] dW_Mtx       = null;	// Estimated Unmixing Matrix.
	private double [][] dS_Mtx       = null;        // Source Matrix.
        private double []   dY_Vec       = null;
	/**/
	//static Connection connection;
	/*********************************************************/
	public static void main (String  argv [])
	{

            if (argv.length == 0){
                System.out.println ("usage: java IcaPreprocess [file] [Component Number] [MaxIteration]");
                return;
            }

            System.out.println (argv [0].toString ());
            iNComp     = Integer.valueOf (argv [1]).intValue ();	// Component nuber to extract
            if (iNComp <= 0) {
                System.out.println ("Error: improper number of components requested");
                return;
            } 
            int iMaxIt = Integer.valueOf (argv [2]).intValue ();    // Max number of iterations
            if (iMaxIt <= 0) {
                System.out.println ("Error: improper number of iterations requested");
                return;
            }

            Preprocess inputIca = new Preprocess (argv [0], iNComp, iMaxIt);
            //new IcaPreprocess ();
	}
        /***********************************************************/
        public Preprocess (int iNCompIn, int iMaxIterIn)
        {
            iNComp     = iNCompIn;
            iMaxIter   = iMaxIterIn;
            //dDataTrain = DinputMtx;            
        }
	/*********************************************************/
	/* Input Data from a text file                           */
	/*********************************************************/
	public Preprocess (String sFileName, int iNCompIn, int iMaxIterIn)
	{
                                                       /*readData gets the X matrix, containing the xi noisy vectors */
            iNComp     = iNCompIn;
            iMaxIter   = iMaxIterIn;
            dDataTrain = readData (sFileName );
            System.out.println ("Readed Data");
            Miscellanea.printMatrix (dDataTrain);

            System.out.println ("Preprocess data.");
            dX_Mtx = PCAPreProcess (dDataTrain);

            dS_Mtx = ICAPreProcess (dX_Mtx);
            Miscellanea.printMatrix (dW_Mtx);

            try {
                DataOutputStream output = new DataOutputStream(new FileOutputStream("W.txt"));
                writeMatrix (output, dW_Mtx);
                output.close();
                output = new DataOutputStream(new FileOutputStream("S.txt"));
                writeMatrix (output, dS_Mtx);
                output.close();
            } catch (Exception e) {
                System.err.println ("Preproces: " +e);
            }	
	}        
	/***********************************************************************/
	/* Input Data from a Double Matrix, presumably comming from databases. */
	/* 0j0 Double != double!!!                                             */
	/***********************************************************************/
	public Preprocess(double[][] DinputMtx, int iNCompIn, int iMaxIterIn)
	{
            iNComp     = iNCompIn;
            iMaxIter   = iMaxIterIn;
            dDataTrain = DinputMtx;

            System.out.println ("Preprocess data.\n" +" iRow:"+dDataTrain.length + " iCol: "+dDataTrain[0].length);

            dX_Mtx = PCAPreProcess (dDataTrain); //DinputMtx);

            dS_Mtx = ICAPreProcess (dX_Mtx);

	}      
	/*********************************************************/
	double [][] readData(String filename)
	{
            int N = 0,
                d = 0;
            BufferedReader in   = null;
            FileReader     file = null;

            try { 
                    file = new FileReader(filename);
                    in   = new BufferedReader(file);
            } catch (Exception e) 
            {
                    System.err.println("Error 0: " +e);
                    System.exit(0);
            }

            try {
                    String sLine = in.readLine ();
                    for(int k = 0; k < sLine.length (); k++) {
                            if(sLine.charAt(k) != ' ') {
                                    iColumn ++;
                                    for (;k < sLine.length (); k++)
                                            if (sLine.charAt (k) == ' ') 
                                                    break;
                                    }
                    }
                    System.out.println ("Cols: " + iColumn);

                    iRow = 1;
                    while (in.readLine () != null) 
                            iRow++;

                    System.out.println ("Rows: " + iRow);

                    in.close ();
            } catch (Exception e) 
            {
                    System.err.println ("Error 1: "+ e);
                    System.exit (0);
            }

            System.out.println ("Reading " + " iRow "+ iRow + " Columnas " + iColumn );

            double [][] dDataTrain = new double [iRow] [iColumn]; 
            try {
                    file = new FileReader(filename);
                    in   = new BufferedReader(file);
                    ParseData(in, dDataTrain);
                    in.close();
            } catch (Exception e) 
            {
                    System.err.println ("Error 2: "+ e);
                    System.exit (0);
            }

            return dDataTrain;
	}
	/*********************************************************************************************/
	/* ParseData: parse a numeric text file that represents a matrix devided in space separated. */
	/* Note: this method works for two columns only.                                             */
	/*********************************************************************************************/
	void ParseData (BufferedReader in, double [][] Data) throws Exception
	{
            int i,
                j,
                k;
            String sLine;
            StringBuffer word = new StringBuffer();

            for(i=0; i< iRow; i++) {
                sLine = in.readLine(); 
                k = 0;
                //System.out.println ("sLine "+ sLine);
                for(j=0; j < iColumn; j++) {
                    word.setLength(0);

                    for(; k < sLine.length (); k++)
                        if(sLine.charAt(k) != ' ') 
                            break;

                    for(; k < sLine.length (); k++) {
                        char c = sLine.charAt(k);
                        if(c == ' ') 
                            break;
                        if (c == '\t')
                            break;
                        word.append (sLine.charAt(k));
                    }

                    Data [i][j] = Double.valueOf(word.toString()).doubleValue();
                    //System.out.print (Data [i ][ j]+"\t");
                }
                //System.out.println ();
            }
        } /*end ParseData method */
        /********************************************************************************************/
	public double [][] PCAPreProcess (double [][] dDataTrain)
	{   
            double [][] dWhiteOut = null;

            System.out.println ("center data");
            dCentMtx   = getCenteredData (dDataTrain);
            
            System.out.println ("scale data");
            dScaledMtx = getScaledData   (dCentMtx);
            
            System.out.println ("white data");
            SingularValueDecomposition mtxE = getEigenValues (dScaledMtx, iNComp);
            
            Matrix mtxEigenValues = mtxE.getS ();
            System.out.println ("EigenValues");
            //mtxEigenValues.print (10, 4);
            dEigenValues = mtxEigenValues.getArrayCopy ();
//            int iRow = mtxEigenValues.getColumnDimension();
//            dEigenValues     = new double [iRow][iRow];
//            for (int i = 0; i < iRow; i++) {
//                dEigenValues [i][i] = mtxEigenValues.get(i, i);
//            }

            Matrix mtxDVals2 = getEigen2Values (mtxEigenValues);
            System.out.println ("EigenVectors");
            //mtxDVals2.print (10, 4);
            
            dKMtx            = getKMatrix (mtxE, mtxDVals2, iNComp);
            
            dWhiteOut        = getWhiteData    (dScaledMtx, dKMtx, iNComp);

            return dWhiteOut;

	}/* end PCAPreProcess */ 
        /********************************************************************************************/
        public double [][] ICAPreProcess (double [][] dX_MtxIn)
        {
            double [][] dS_MtxOut = null;
            
            System.out.println ("fastICA. Type: deflation. Negentrphy finction: EXP.");
            dW_Mtx       = getExpDefICA (dX_MtxIn);

            Matrix mtx_W = new Matrix (dW_Mtx);
            dWMtx = mtx_W.getArrayCopy ();
            Matrix mtx_X = new Matrix (dX_Mtx, iRow, iNComp);
//		//mtx_X.print (10, 4);
            Matrix mtx_S = mtx_W.times (mtx_X.transpose ());
            mtx_S = mtx_S.transpose ();
//		//mtx_S.print (10, 4);
            dS_MtxOut  = mtx_S. getArray (); 
            return dS_MtxOut;
        }
        /********************************************************************************************/
	public void writeMatrix (DataOutputStream outputParam, double [][] DataIn)
	{
            int i,
                j;

            int iRowNum = DataIn .length ;
            int iColNum = DataIn [0].length;
            try {	
                for (i = 0; i < iRowNum; i++) {
                    for (j = 0; j < iColNum; j++) {
                        //System.out.print ("["+1+"]["+j+"] "+DataIn [i][j]+"\t");
                        Double dTmp = new Double (DataIn [i][j]);
                        outputParam.writeBytes (dTmp.toString ()+" ");
                    }
                    outputParam.writeBytes ("\n");
                }
            } catch (IOException e) {}
	
	} /* end writeMatrix */       
        /********************************************************************************************/ 
	Matrix setDiag (int iColNum, double [] DataIn)
	{
            Matrix mtxDiag = new Matrix (iColNum, iColNum);
            double [][] ptrDiag = mtxDiag.getArray ();
            //System.out.println ("Columnas  de la Diagonal "+ iColNum);
            for (int i = 0; i < iColNum; i++)
                    ptrDiag [i][i] = DataIn [i];

            return mtxDiag;
	} /* end SetDiag */
        /***********************************************************************************************/
        /* The param row.norm of the R fastICA code is directly involved with the data centering.      */
        /* row.norm: a logical value indicating whether rows of the data matrix                        */
        /*          `X' should be standardized beforehand.                                             */
        /* On R, this function actually centers and scales a vector                                    */ 
        /*                                                                                             */
        /* Checked against R algorithm for option row.norm = TRUE                                      */
        /***********************************************************************************************/
	double [][] getCenteredData (double [][] DataIn)
	{
            int iRowNum = DataIn.length;
            int iColNum = DataIn [0].length;
            double [][] dMtxOut  = new double [iRowNum][iColNum];

            //System.out.println ("mean");
            double [] dTempMean = Miscellanea.getMeanVector (DataIn); 

            //System.out.println ("after mean");
            for (int i = 0; i < iRowNum; i++) {
                for (int j = 0; j < iColNum; j++) {			
                    dMtxOut [i][j] = (DataIn [i][j] - dTempMean [j]);

                    //System.out.print (dMtxOut [i][j]+"\t");
                }
                //System.out.println ();
            }

            return dMtxOut;
	} /* end getCenteredData */
        /***********************************************************************************************/
	double [][] getScaledData (double [][] DataIn) //int iRowNum, int iColNum, double [][] DataIn)
	{
            int iRowNum            = DataIn.length;
            int iColNum            = DataIn [0].length;            
            double [][] dTemp      = Miscellanea.getTransposeMtx (DataIn);
            double [][] dScaleMtx  = new double [iColNum][iRowNum]; 

            System.out.println ("Scaled: CenterData");
            dScaleMtx = getCenteredData (dTemp);
            System.out.println ("RootMeanVect");
            double [] dTempRootM = Miscellanea.getRootMeanVect (iColNum, iRowNum, dScaleMtx); 

            System.out.println ("\nScaled: Scaling");
            for (int i = 0; i < iColNum; i++) {
                for (int j = 0; j < iRowNum; j++) {			
                    dTemp [i][j] = (dScaleMtx [i][j]/dTempRootM [j]); //<<<<<<<<<<<<<<<<<<<<<<0j0

                    //System.out.print (dTemp [i][j]+"\t");
                }
                //System.out.println ();
            }

            return Miscellanea.getTransposeMtx (dTemp);
	} /* end getScaledData */
        /********************************************************************************************/
        SingularValueDecomposition getEigenValues (double [][] dMtxIn, int iNCompIn)
        {
            int iRowNum = dMtxIn.length;
            int iColNum = dMtxIn [0].length;
                        
            Matrix mCentered = new Matrix (Miscellanea.getTransposeMtx (dMtxIn));
            //mCentered.print (10, 4);

            System.out.println ("Covar mtx");
            Matrix mtxV = mCentered.times (mCentered.transpose ());
            double [][] ptrV = mtxV.getArray ();
            for (int i = 0; i < iColNum; i++) {
                for (int j = 0; j < iColNum; j++) {			
                    ptrV [i][j] = (ptrV [i][j]/iRowNum);       //<<<<<<<<<<<<<<<<<<<<<<0j0

                    //System.out.print (ptrV [i][j]+"\t");
                }
                //System.out.println ();
            }
            //aca mtxV.print (10, 4);

            SingularValueDecomposition mtxE = mtxV.svd ();
            
            return mtxE;
        }
        /********************************************************************************************/
        Matrix getEigen2Values (Matrix mtxEigenValues)
        {
            double ptrEigenValues [][] = mtxEigenValues.getArray ();
            int iCol = ptrEigenValues [0].length;
            
            for (int i = 0; i < iCol; i++) {
                double dTmp = Math.sqrt (Math.abs (ptrEigenValues [i][i]));
                ptrEigenValues [i][i] = (dTmp != 0)? 1/dTmp: 0;                   //<<<<<<<<<<<<<<<<<<0j0
            } 
            
            return new Matrix (ptrEigenValues);
        }
        /********************************************************************************************/
        double [][] getKMatrix (SingularValueDecomposition mtxE, Matrix mtxDVals, int iNCompIn)
        {
            System.out.println ("Mtx K");
            Matrix mtxEU    = mtxE.getU ();
            Matrix mtxKTemp = mtxDVals. times (mtxEU.transpose ());
            int iColDim     = mtxKTemp.getColumnDimension();
            //mtxKTemp.print (iColDim, 4);

            double [][] ptrKTmp = mtxKTemp.getArrayCopy ();
            double [][] ptrK    = ptrKTmp;
            
            System.out.println ("Adjusted K; nCompIn:" +iNCompIn + " iColDim: " + iColDim);
            for (int i = 0; i < iColDim; i++) {
                for (int j = 0; j < iColDim; j++) {
                    //if (i >= (iColDim - iNCompIn)) 
                    if ( i < iNCompIn)	//<<<<<<<<<<<<<<<0j0	
                            ptrK [i][j] = -ptrK [i][j];
                    else
                            ptrK [i][j] = 0;

                    //System.out.print (ptrK [i][j]+"\t");
                }
                //System.out.println ();
            }

            return ptrK;
        }
        /********************************************************************************************/
        /* getWhiteData: Whitening procedure of a centered vector, as descrived in doc p. 15        */
        /* input params:                                                                            */
        /********************************************************************************************/
        public double [][] getWhiteData (double [][] dMtxIn, double [][] dKMtx, int iNCompIn) //throws Exception
        {
            int iRow = dMtxIn.length;
            int iCol = dMtxIn [0].length;
                        
            Matrix      mtxIn  = new Matrix (Miscellanea.getTransposeMtx (dMtxIn));
            System.out.println ("KRow:" + dKMtx.length+" KCol: " +dKMtx[0].length+" Dif: "+(iCol -iNCompIn));
            Matrix      mtxK   = new Matrix (dKMtx);
            Matrix      mtxPCA = mtxK.times (mtxIn);
            mtxPCA             = mtxPCA.transpose ();
            double [][] dTmpMtx = mtxPCA.getArray ();
            double [][] dPCAMtx = new double [iRow][iNCompIn];
            
            for (int i = 0; i < iRow; i++) {
                for (int j = 0; j < iNCompIn; j++) {
                    dPCAMtx [i][j] = dTmpMtx [i][j];
                    //System.out.print (dPCAMtx [i][k++]+" ");
                }
                //System.out.println ();
            }
            
            return dPCAMtx;
        }             
        /***********************************************************************************************/
        /* getWeigthMtx: get a random weight matrix in order for ICA to start to iterate.              */
        /*                                                                                             */
        /* Note: This ougth to be done with a neuron in the future.                                    */
        /***********************************************************************************************/
	Matrix getWeigthMtx (int iNCompIn, String sStage, Matrix mtxWC)
	{	
            Matrix mtxW = null;

            if (sStage == "I" || sStage == "i") {  // I stands for initial stage.
                mtxW = new Matrix (iNCompIn, iNCompIn); 
                mtxW = mtxW.random (iNCompIn, iNCompIn);
                //mtxW.print (10, 4);
            } 
            if (sStage == "C" || sStage == "c") { // C stands for Continuation stage
                if (mtxWC == null) {
                    System.out.println ("Error: Weigth matrix pointer = null");
                    return null;
                }
                mtxW = mtxWC;
            }	

            SingularValueDecomposition mtx_sW = mtxW.svd ();
            Matrix mtx_dW = mtx_sW.getS ();
            double [][] ptr_dW = mtx_dW.getArray ();
            //mtx_dW.print (10, 4);
            for (int i = 0; i < iNCompIn; i++) {
                if (ptr_dW [i][i] == 0) {
                    ptr_dW [i][i] = 0; 
                    System.err.println ("Warning: Attempt to divide by zero.");       //<<<<<<<<<<<<<<0j0
                    continue;
                }
                ptr_dW [i][i] = 1/ptr_dW [i][i]; 
            }
            //mtx_dW.print (10, 4);
            Matrix mtx_uW   = mtx_sW.getU ();
            Matrix mtx_Temp = new Matrix (iNCompIn, iNCompIn);
            mtx_Temp        = mtx_uW.times (mtx_dW);
            mtx_Temp        = mtx_Temp.times (mtx_uW.transpose ());
            mtxW            = mtx_Temp.times (mtxW);
            //aca mtxW.print (10, 4);

            return mtxW;
	} /* end getWeigthMtx */
        /***********************************************************************************************/
        /* Estimates an unmixing matrix W given the negentrophy function W.                            */
        /***********************************************************************************************/
	Matrix expNegEntrophy (int iRowNum, int iNCompIn, Matrix mtx_wx, String vx)
	{
            Matrix mtxTemp = new Matrix ( mtx_wx.getColumnDimension (), mtx_wx.getRowDimension ());
            mtxTemp = mtx_wx;

            double [][] ptr_Temp  = mtxTemp.getArray ();
            double [][] ptr_wx    = mtx_wx.getArray ();
            double [][] d_expwx   = new double [iNCompIn][iRowNum];

            for (int i = 0; i < iNCompIn; i++) {
                for (int j = 0; j < iRowNum; j++) {
                    Math.pow (ptr_Temp [i][j], 2);
                    if (vx == "v1" || vx == "V1")
                        d_expwx [i][j] = ptr_wx [i][j]  * Math.exp (-(1/2)*ptr_Temp [i][j]);
                    else if (vx == "v2" || vx == "v2")
                        d_expwx [i][j] = (1 - ptr_Temp [i][j]) * Math.exp (-(1/2)*ptr_Temp [i][j]);
                    //System.out.print (ptr_Temp [i][j]+" ");			
                }
                //System.out.println ();
            }
            //mtx_wx.print (10, 4);
            Matrix mtx_expwx = Matrix.constructWithCopy (d_expwx);

            return mtx_expwx;
	} /* end expNegEntrophy */
        /***********************************************************************************************/
        /* getExpDefICA: ICA agorithm with the following features:                                     */
        /* negentrophy function type: Exponential.                                                     */
        /* Component extraction:      Deflation i.e. get one ICA componet at the time.                 */
        /***********************************************************************************************/
	double [][] getExpDefICA (double [][] DataIn)
	{
            int    i,
                   j;
            int iRowNum = DataIn.length;
            int iColNum = DataIn [0].length;

            if (iNComp > iColNum) {
                System.out.println ("Error: Number of columns bigger than number of components");
                return null;
            }

            System.out.println ("Initial matrix W:");
            Matrix mtxW         = getWeigthMtx (iNComp, "I", null);
            Matrix mtxX         = new Matrix   (DataIn, iRowNum, iNComp);
            Matrix mtx_XTemp    = new Matrix   (DataIn, iRowNum, iNComp);
            mtx_XTemp = mtx_XTemp.times (1/mtx_XTemp.getRowDimension ());

            Matrix mtxW1 = null;

            for (int it = 0; it < iMaxIter; it++){
                //mtx_wx.print (10, 4);
                Matrix mtx_wx   = mtxW.times      (mtxX.transpose ());
                Matrix mtx_gwx  = expNegEntrophy  (iRowNum, iNComp, mtx_wx, "v1");
                Matrix mtxV1    = mtx_gwx.times (mtx_XTemp);

                Matrix mtx_g_wx  = expNegEntrophy  (iRowNum, iNComp, mtx_wx, "v2");
                Matrix mtx_gTmp  = mtx_g_wx.transpose ();
                double [] dgMean = Miscellanea.getMeanVector (mtx_gTmp.getArray ());
                Matrix mtxV2     = setDiag (mtx_gTmp.getColumnDimension (), dgMean);
                mtxV2 = mtxV2.times (mtxW);

                mtxW1 = mtxV1.minus (mtxV2);
                mtxW1 = getWeigthMtx (iNComp, "C", mtxW1);
                System.out.print ("+");
                //aca mtxW1.print (10, 4);
                //mtx_wx.print (10, 4);
                //mtx_XTemp.print (10, 12);
            }
            System.out.println ("+"); //"W1");
            //aca mtxW1.print (10, 4);

            return mtxW1.getArray ();
	} /* end getScaledData */
        /*********************************************************/
	public double [][] getSMatrix ()
	{
            if (dS_Mtx == null)
                return null;

            return dS_Mtx;
	} /* end getSMatrix */   
        /*********************************************************/
        public double [] getYVect ()
        {
            return dY_Vec;
        }
        /*********************************************************/
        public double [][] getIndependentComponents (double [][] dInputMtx)
        {
            Matrix mtx_S = new Matrix (dInputMtx);
            Matrix mtx_W = new Matrix (dW_Mtx);
            Matrix mtx_Y = mtx_W.times(mtx_S.transpose());
            mtx_Y = mtx_Y.transpose();

            double [][] dYMtx = mtx_Y.getArrayCopy();
            
            return dYMtx;
        }
        /************************************************************/
        public double [][] getEigenValues ()
        {
            return dEigenValues;
        }
        /*************************************************************/
        public double [][] getKMatrix ()
        {
            return dKMtx;
        }
        /*************************************************************/
        public double [][] getWMatrix ()
        {
            return dWMtx;
        }
} /* end IcaPreprocess class */








